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EXINGTON 


MASSACHUSETTS 


ABSTRACT 


In  this  pahcr,  we  develop  a  set  of  conditions  under  which  a  sequence 
is  uniquely  specified  by  the  phase  or  samples  of  the  phase  of  its 
Fourier  transform,  and  a  similar  set  of  conditions  under  which  a  sequence 
is  uniquely  specified  by  the  magnitude  of  its  Fourier  transform.  These 
conditions  are  distinctly  different  from  the  maximum  or  minimum  phase 
conditions,  and  are  applicable  to  both  one-dimensional  and  multi-dimensional 
sequences.  Under  the  specified  conditions,  we  also  develop  several 
algorithms  which  may  be  used  to  reconstruct  a  sequence  from  its  phase  or 
magnitude.  As  a  potential  application  area,  the  results  of  this  i-aj)er 
are  applied  to  the  blind  deconvolution  problem  of  digital  images  blurred 
by  a  symmetric  point  spread  function. 
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I.  INTRODUCTION 

For  both  continuous-time  and  discrete-time  signals,  the  magnitude 
and  phase  of  the  l-ourier  transform  are,  in  general,  independent  functions, 
i.e.,  the  signal  cannot  be  recovered  from  knowledge  of  either  alone. 

Under  certain  conditions,  however,  relationships  exist  between  these 
components.  For  example,  when  the  signal  is  a  minimum  phase  or  maximum 
phase  signal  both  the  log  magnitude  and  phase  can  be  obtained  from  the 
other  through  the  Hilbert  transform.  This  relationship  has  been  exploited 
in  a  variety  of  ways  in  many  fields  including  network  theory,  communications 
and  signal  processing  (l,2,3j. 

In  this  paper  we  develop  a  set  of  conditions  under  which  a  discrete¬ 
time  sequence  is  completely  specified  to  within  a  scale  factor  by  the 
phase  of  its  Fourier  transform,  without  the  restriction  of  minimum  or 
maximum  phase,  and  propose  several  algorithms  for  implementing  the 
reconstruction  of  a  signal  from  the  phase  of  its  Fourier  transform.  In 
Section  II  we  consider  the  case  in  which  the  phase  is  specified  at  all 
frequencies,  and  in  Section  in  the  case  where  the  phase  is  specified  at 
a  discrete  set  of  frequencies.  Algorithms  for  implementing  the  reconstruction 
are  developed  in  Section  V.  In  Section  IV,  we  develop  a  different  set 
of  conditions,  again  without  the  restriction  of  minimum  or  maximum 
phase,  in  which  a  discrete-time  sequence  is  completely  specified  by  the 
magnitude  of  its  Fourier  transform.  In  Section  VI,  we  extend  the  results 
of  Sections  II,  III,  IV  and  V  to  the  case  of  multi-dimensional  sequences 
and  illustrate  an  application  of  the  results  of  this  paper  to  blind 


deconvolution. 


II.  UNIQUUNLSS  OF  A  SLQUENCF  WITH  A  PllASF  FUNCTION  SPFCIFILU  AT 


ALL  FREQULNCILS 


The  sequences  that  we  consider  are  real  with  rational  z-transf orms . 

Since  we  are  interested  in  conditions  under  which  the  sequence  can  be 

uniquely  specified  by  the  phase  of  its  Fourier  transform,  the  lourier 
transform  is  assumed  to  converge  i.e.,  the  region  of  coiv'crgence  of  the 
z-transform  includes  the  unit  circle. 

For  such  sequences,  we  first  show  that  a  finite  length  sequence  is 
uniquely  specified  by  the  phase  of  its  Fourier  transform  if  its  z- 
transform  has  no  zeroes  in  reciprocal  pairs  or  on  the  unit  circle.^  More 

specifically,  denoting  the  phase  of  x[n]  and  y[n]  by  and 

0y(u))  respectively,  we  demonstrate  the  following: 


Theorem  1:  Let  x[n]  and  y[n]  be  two  finite  length  sequences 

whose  z-transforms  have  no  zeroes  in  reciprocal 

pairs  or  on  the  unit  circle.  If  0  (w)  =  O(i^) 

X  y 

for  all  u,  then  x [n]  =  Ay[n]  for  some  positive 

constant  A.  If  tan  0  (w)  =  tan  0  (u)  for  all  co, 

X  y 

then  x[nj  =  Ay[n]  for  some  real  constant  A. 


To  demonstrate  the  validity  of  Theorem  1,  we  note  first  of  ali, 
that  if  a  finite- length  sequence  g[n]  with  z-transform  G(z)  has  a  phase 
which  is  zero  or  tt  for  all  w,  then  g[n]  is  an  even  sequence,  i.e.,  g[n]= 
g[-n],  and  consequently,  if  G(z)  has  a  zero  at  z=z^,  then  G(z)  must  also  have 


Since  we  are  considering  only  sequences  which  are  real,  zeroes  occur  in 
complex  conjugate  pairs.  In  the  following  discussions,  this  symmetry  is 
implicitly  assumed,  particularly  in  reference  to  reciprocal  zeroes. 
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I  *n 


a  zero  at  2=1/2^.  Now  assume  that  x(n]  and  y[n]  both  satisfy  the  conditions 
of  Theorem  1  and  define  g[n]  as 


so  that 


g[n]=x[n]*y[-n] 

GC2)=X(z)Y(1/z) 


If  or  if  tan  0^(0)),  then  the  phase  of  g[n] 

is  zero  or  tt.  Therefore,  g[n]  is  an  even  sequence.  Since  the  zeroes 

of  G(2}  are  collectively  the  zeroes  of  X(z)  and  Y(l/z),  if  X(Zq)=0 

then  either  X(1/Zq)=0  or  Y(Zq)=0.  However,  because  the  conditions  of 

Theorem  1  exclude  reciprocal  zeroes  or  zeroes  on  the  unit  circle,  X(Zq) 

and  X(1/Zq)  cannot  both  be  zero.  Thus,  if  X(2j^)=0  then  Y(2q)=0  and 

vice  versa,  i.e.,  the  zeroes  of  X(z)  and  Y(z)  are  identical.  Consequently, 

2 

since  g[n]  is  an  all  zero  sequence. 


X(2)=AY{2)  (3)  a 

or  x[n]=Ay[n]  (3)b 

Combining  equations  1  and  3b,  we  have 

g[n]=Ax[n] *x[-n]  (4) 

Since  the  phase  of  x[n]*x[-n]  is  always  zero,  if  0  (co)  =  0  Cw)  then  the 

x  y 

phase  of  g[n]  is  zero  so  A  must  be  a  positive  constant.  If  tan  = 

tan  0y(co)  then  the  phase  of  g[n]  is  zero  or  tt  so  A  must  be  real. 

An  interpretation  of  Theorem  1  is  suggested  by  the  observation  that 

for  a  rational  z-transform,  in  general,  a  zero  at  z=Zq  and  a  pole  at 
★ 

z=1/Zq  contribute  the  same  phase  but  different  magnitude  to  the  Fourier 


When  we  refer  to  a  sequence  as  an  all-zero  (all-pole)  sequence  this 
should  be  interpreted  to  mean  that  the  z-transform  has  only  zeroes 
(poles)  except  possibly  at  z=0  or  z='”. 
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transform.  I'hus  with  phase  information  alone,  there  is  an  inherent  ambigui  y 

in  the  t-transform  in  the  sense  that  a  zero  (polcj  at  ---jj  associated  with  the 

original  seciuence  can  potentially  only  be  identified  from  the  given  phase 

★ 

as  either  a  zero  (.polej  at  or  a  pole  (zero)  at  z=l/z^  and  this 

ambiguity  cannot  be  further  resolved  without  additional  information  or 
conditions.  The  finite  length  condition  in  Theorem  1  resolves  this  ambiguity 
by  restricting  the  z-transform  to  have  only  zeroes  except  ])ossibly  at  z=0 
or  at  z=°°.  The  additional  condition  that  the  z-transform  have  no  zeroes 
in  reciprocal  pairs  eliminates  the  possibility  of  zero  phase  components 
in  the  z-transform  which,  of  course,  could  never  be  recovered  from  phase 
information  alone.  The  conditions  in  Theorem  1  also  eliminate  the 
possibility  of  zeroes  on  the  unit  circle.  While  the  Theorem  can  be 
modified  to  allow  for  the  possibility  of  zeroes  on  the  unit  circle,  the 
result  becomes  somewhat  more  complicated  and  we  have  chosen  not  to 
include  this  additional  generality. 

Although  Theorem  1  requires  that  x[n]  be  an  all-zero  sequence,  a 
dual  to  Theorem  1  can  be  formulated  for  an  all-pole  sequence.  Specifically, 
let  x[n]  denote  the  convolutional  inverse  of  a  sequence  x[n],  i.e. 

x[n]  *  x[n]  =  -S[nl  (5) 

Then : 

Theorem  2:  Let  x[n]  and  y[nl  be  two  sequences  whose  z-transforms 

have  no  poles  in  rcci()rocal  pairs,  and  which 

have  finite  duration  convolutional  inverses.  If  0  (co)  = 

0^  (Lo)  for  all  w,  then  x|n]  =  Ay[n]  for  some  positive 

constant  A.  If  tan  0  (la)  =  tan  0  (la)  for  all  u),  then 

X  y 

x[nJ=Ay[n|  for  some  real  constant  A. 


Theorem  2  follows  directly  from  Theorem  1.  Since  the  phase  of  the 
Tourier  transform  of  x[n]  is  specified  by  the  j^hase  of  the  Fourier 
transform  of  x[n],  x[n|  is  uniquely  specified  to  within  a  positive 
scale  factor  by  the  phase  of  the  Fourier  transform  of  x[n],  by  virtue 
of  Theorem  1.  Then,  x[n)  is  uniquely  determined  from  the  reciprocal 
of  the  z-transform  of  x[n]. 

In  Section  IV  we  will  consider  a  number  of  numerical  algorithms 
which  can  be  implemented  on  a  digital  computer  for  reconstructing  a 
sequence  from  its  phase  under  the  condition  of  Theorem  1  or  Theorem  2. 

At  this  point,  however,  we  discuss  a  conceptual  algorithm  which  may 
potentially  have  a  practical  implementation  but  which,  more  importantly,  serves 
to  lend  insight  into  Tiieorem  1  and  2.  We  outline  the  algorithm  under 
the  conditions  of  Theorem  1  since  it  is  easily  modified  for  the  conditions 
of  Theorem  2. 

Let  0^(w)  denote  the  specified  phase  function  from  which  the  sequence 
is  to  be  reconstructed  and  the  associated  unwrapped  phase  (5).  From 

the  conditions  of  Theorem  1,  X(z)  is  restricted  to  be  of  the  form 

^0  ^'l  -1  ^2 

X(z)=Cz  (l-a.z  ')  n'"  (1-b,  z)  (0) 

k=l  k=l 

'1,  \\ 

k  and 
Step  1: 

The  algebraic  sign  of  C  is  obtained  from  using  the  fact  that 

0^(0}  is  zero  if  and  only  if  C  is  positive  (3) .  The  value  of  n^  in 


with  C  real,  n^  an  integer,  |aj^| 


★ 

<1  for  all  k  and  aj^^  b^  for  any 


eq.  (6)  is  obtained  from  the  unwrapped  jjhasc  as 


n 

o 


i  -  o/»)l 


Step  2: 

I'roin  the  unwrapped  phase  function  and  the  value  of  n^  obtained  in 
Step  1,  a  new  pliase  function  is  specified  as 


n  cj  -  0  (0) 
o  ^ 


(8) 


Using  the  Hilbert  transform,  a  minimum  phase  sequence  can  be 

specified  which  has  the  phase  0^(w).  The  z-transform  (z)  of  fn] 


is  given  by  (3,4) 


min 


(t)  = 


*'l  - 1 

n  ) 


2  *  '1 

r  (i-b.  z  ■") 

k=l 


(9) 


where  the  coefficients  and  bj^  are  identical  to  those  in  eq.  (6). 

Since  pole-zero  cancellations  cannot  occur  in  eq.  (9)  by  virtue  of  the 

condition  in  Theorem  1  which  implies  that  for  any  k  or  Z,  the  coefficients 

in  eq.  (6)  can  be  obtained  from  the  zeroes  of  X  .  (z)  and  the  coefficients 

^  ^  min 

* 

b,  ,  and  thus  b,  in  eq.  (6)  can  be  obtained  from  the  poles  of  X  .  (z) . 

The  condition  in  Theorem  1  that  there  are  no  zeroes  in 
reciprocal  pairs  ensures  that  there  are  no  pole-zero  cancellations  in  eq. 

(9).  If  the  original  sequence  has  reciprocal  zeroes,  then  the  algorithm 
above  may  still  be  applied  to  recover  all  but  those  zeroes  in  X(z) 
which  are  in  reciprocal  pairs. 
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1>1S0KI  II-  I  Ul  ijlll-MIll.S 

In  Tlicorfiii.s  1  and  J  we  assumed  that  the  phase  function  was  specified 
at  all  frequencies.  A  similar  set  of  Theorems  can  be  stated  if  the  phase 
is  specified  at  a  sufficient  number  of  discrete  frequencies.  As  in 
Section  11  wo  assume  the  sequences  are  real  with  rational  z-transforms 
with  a  region  of  convergence  that  includes  the  unit  circle. 

Then : 

Theorem  5:  Let  x[n]  and  y[n]  be  two  finite  length  sequences 

which  are  zero  outside  the  interval  0<n<.\-l  with 

z-transforms  which  ha\ e  no  zeroes  in  reciprocal 

pairs  or  on  the  unit  circle.  If  G  (u)  =  0  (w) 

X  y 

at  (N-1)  distinct  frequencies  greater  than  zero 
and  less  than  "n.then  xin)  =  Ay[Ti]  for  some 
positive  constant  A.  If  tan  0^(cjJ  =  tan  (‘^') 
at  (N-1)  distinct  frequencies  greater  than  zero 
and  less  than  ti,  then  x[n]  =  Ay[nJ  for  ^■'■'mc  real 
constant  A. 

The  validity  of  Tlteorem  5  follows  in  an  almost  identical  manner 
to  that  of  Theorem  1.  Specifically,  consider  a  finite  length  sequence 
g[nj  for  which  g[n]=0  outside  the  interval  (-iN+ 1 ) <n<  (N- 1 )  .  Let 

I  1 

G(w)  =  I  G  (w)  I  e  ^  denote  the  Fourier  transform  of  g[nj  with 
zero  or  tt  at  N-1  discrete  frequencies  Wj , w,, , .  .  .  ^  between  o'=0  and 

u=Tr ,  i  .  e .  , 

'^More  generally,  x[n]  need  only  he  zero  outside  any  finite  interval 
of  length  N.  This  added  generality,  however,  is  not  considered  in 
order  to  simplify  the  following  discussions. 
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(lOJa 


0  (u),  )  =  0  or  Ti  for 
g  k 

k-1. 

with 

^  for 

and 

0<Lo  <1!  for  k=l,... 

K 

,.N-1 

Then , 

t; ((jj|^)  is  real  and 

N-1 

lm[t:(tj,)J=  L  g[nj 

n=-N+l 

sin 

Or,  equivalently 


N'-l 

Z 

n=  1 


sin  nwj^=0 


k=l,2. . .N-1 

0<W|,<TI 


liquation  (12)  imjilies  that  (5) 

[g[n]-g[-n]]  =0  n=l,...N-l 


(.0)h 

(10)c 


(11) 


(12) 


(13) 


i.e.,  g[n]  is  an  even  sequence.  Now,  consider  two  sequences  x[n]  and 
y[n]  satisfying  the  conditions  of  Theorem  3  and  having  the  same  phase  at 
(N-1)  distinct  frequencies  between  zero  and  tt.  As  with  Theorem  1,  we 
form  the  sequence 

g[n]=x[n]*y[-n]  (14) 

Since  g[n]=0  outside  the  interval  (-N+l)<n< (N-1)  and  0  (w.)  satisfies 
eqs.  (10),  then  g[nj  is  an  even  sequence.  For  reasons  identical  to 
those  used  in  justifying  Theorem  1,  it  then  follows  that 

x[n]=Ay[nJ  (15) 

where  A  is  a  pf^sitive  constant  if  the  phase  samples  of  x[n]  and  y[n]  ar' 
equal  and  a  real  constant  if  tangents  of  the  phase  samples  are  equal. 

Although  Theorem  3  requires  that  x[n]  be  an  all-zero  sequence,  a 
dual  to  Theorem  3  for  an  all-pole  sequence  is  easily  formulated  in  terms 
of  the  convolutional  inverse.  Specifically; 
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Theorem  4 : 


Let  x[n]  and  /[n]  be  two  sequences  whose  z-transforms 

have  no  poles  in  reciprocal  pairs,  and  which 

have  convolutional  inverses  that  are  zero  outside  the 


interval  0<n<N-l.  If  0  (wj  =  0  (cu)  at  CN-1)  distinct 
-  -  X  y 

frequencies  greater  than  zero  and  less  than  tt  then  x[n]  = 
Ay[n]  for  some  positive  constant  A.  If  tan  0^(cj)  = 
tan  0^(0))  at  (N-1)  distinct  frequencies  greater  than 
zero  and  less  than  tt,  then  x[n]=Ay[n]  for  some  real 
constant  A. 

IlTeorem  4  follows  from  Theorem  3  in  the  same  manner  that  Theorem  2 
follows  from  Theorem  1„ 

It  should  be  noted  from  Theorem  1,  that  if  the  phase  is  specified 
at  all  frequencies,  then  the  interval  outside  of  which  the  sequence  is 
zero  need  not  be  known  whereas  from  Theorem  3,  if  the  phase  is  specified 

N-1  frequencies  then  this  interval  must  also  be  specified.  In  both 
cases,  except  for  a  positive  scale  factor,  the  phase  in  the  case  of 
Theorem  1  or  phase  samples  in  the  case  of  Theorem  2  uniquely  specify  a 
sequence  within  the  class  of  sequences  which  have  no  zeroes  in  reciprocal 
pairs  or  on  the  unit  circle.  However,  if  the  additional  constraint  that 
x[0]i^0  is  included  in  the  conditions  of  Theorem  5,  then  as  we  now  show, 
the  phase  samples  uniquely  speicify  a  sequence  within  a  broader  class  of 
sequences.  Specifically; 

Theorem  5:  Let  x[nj  be  a  sequence  that  is  zero  outside  the 

interval  0<n<N-l  with  x[0]^0  and  which  has  a 
z-transform  with  no  zeroes  in  reciprocal  pairs 
or  on  tlie  unit  circle.  Let  y[n]  be  any  sequence 
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ivhicti  is  ^ero  outside  the  interval  0<n<N-].  If 
0^  (uij  =0^ (go)  at  (N-IJ  distinct  frequencies  greater 
tlian  zero  and  less  than  ii,  then  y[n]=Ax[nl  for 
some  jiositive  constant  A.  If  tan  0  (go)  =  tan  0  (ojJ 

y  ^ 

at  (N-1)  distinct  frequencies  greater  than  zero  and 
less  than  n,  then  x[n]=Ay[n]  for  some  real  constant  A. 


To  demonstrate  the  validity  of  Theorem  5,  we  first  form  the  sequence 
g[nj =x[nj *y [-n] .  As  discussed  in  Section  III,  since  g[n]=0  outside 
the  interval  (-N+l)  <n<  (N-1)  and  satisfies  equation  (10),  g[n]  is  an 
even  sequence.  Now  let  N^-1  represent  the  location  of  the  last  non-zero 
point  in  x[nj,  i.e.,  x[n]=0  for  n>Nj  and  x[Nj-l]?!0.  Then 

N-1  j 

G(z)  =  X(z)Y(z'^)  =  X  xfnjz'”  Z  y[nlz'"  d^) 

n=()  n=() 

Since  g[n]  is  even  and  x[0]?^0,  y[n]=0  for  n>Nj^  so  that  the  number  of 
zeroes  of  y[n]  is  less  than  or  equal  to  the  number  of  zeroes  of  x[n] . 
Now,  for  reasons  identical  to  those  used  in  justifying  Theorem  1,  if 
g[n]  is  an  even  sequence  and  if  x[n]  has  no  zeroes  in  reciprocal  pairs, 
then  for  each  zero  of  x[n],  y[n]  must  also  have  the  same  zero.  Even 
though  y[n]  is  not  restricted  to  the  class  of  sequences  with  no  zeroes 
in  reciprocal  pairs  or  on  the  unit  circle,  from  our  previous  result, 
y[n]  cannot  have  more  zeroes  than  x[n],  and  therefore,  y[n]=Ax[n]. 

For  reasons  identical  to  those  used  in  justifying  Theorem  1,  A  is  a 
positive  constant  if  the  phase  samples  of  x(n]  and  y[n]  are  equal 
whereas  A  is  a  real  constant  if  the  tangent  of  the  jihase  samples  are 


Although  Theorem  S  requires  that  the  sequence  he  an  all -zero 
sequence,  a  dual  to  Theorem  3  for  all  all-pole  sequence  is  easily 
formulated  in  terms  of  tlie  convolutional  inverse.  Specifically; 


Theorem  b:  Let  x[n]  be  a  sequence  whose  z-transform  has  no 

poles  in  reciprocal  pairs  or  on  the  unit  circle, 
and  whose  convolutional  inverse  is  zero  outside 
the  interval  0<n<iN-l  and  non- zero  at  n=0.  Let 
y(n]  be  any  sequence  whose  convolutional  inverse 
is  zero  outside  the  interval  0<n<N-l.  If  0  (u))  = 

y 

0^(oj)  at  (N-1)  distinct  frequencies  greater  than 

zero  and  less  than  it,  then  y[n]  =  Ax(n]  for  some 

positive  constant  A.  If  tan  0  (a))=tan  0  (w)  at 

y  X 

(N-1)  distinct  frequencies  greater  than  zero  and 
less  than  tt,  then  x[n]=Ay[n]  for  some  real  constant  A. 

Theorem  6  follows  from  Theorem  5  in  the  same  manner  that  Theorem  2 
follows  from  Theorem  1. 

IV.  UNIQUENESS  OF  A  SF.QULNCE  WITH  A  SPECIFIED  MAGNITUDE  FUNCTION 
In  Section  II,  several  sets  of  conditions  are  presented  which 
establish  a  uniqueness  between  a  sequence  and  its  phase  function. 

Unlike  the  case  for  minimum  or  maximum  phase  sequences,  there  is  no 
dual  statement  of  uniqueness  between  a  sequence  and  its  magnitude 
function  under  the  same  set  of  conditions.  However,  under  a  different 
set  of  conditions  a  sequence  is  uniquely  specified  to  within  a  sign  and 
a  time  shift  by  the  magnitude  of  the  Fourier  transform.  The  conditions 
are  embodied  in  the  following  Theorems; 
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riiool'Clu  7;  Let  x[nl  and  )'[ul  I'c  twi)  SL‘i(iiencfs  whose  z-transl'orms 

coutaiii  no  reciprocal  pole-zero  pairs  and 
wiiich  luivc  all  iH>lcs,  not  at  z=-’,  inside  t  lie  unit 
circle  and  all  zeroes,  not  at  z-i),  outside  the  unit 
circle.  11'  the  iiiaj;n  i  tiides  of  the  I'ourier  transforms 
of  \ln|  and  >'ln]  are  eipia  1  then  x|n|  =  +  y|n  +  m]  for 
some  integer  m. 

A  dual  to  this  i'heorem  is: 

Theorem  8:  Lot  x[n]  and  y[n|  he  two  se(|ucnccs  whose  z-t ransforms 

contain  no  reciprocal  pole-zero  pairs  and 
which  have  all  poles,  not  at  z=(),  outside  the  unit 
circle  and  all  zeroes,  not  at  z='",  inside  the  unit 
circle.  If  the  magnitudes  of  the  I'ourier  transforms 
of  x[nl  and  y[nl  are  equal  then  x[nl  =  +  yln+m)  for 
some  integer  ra. 

Since  the  justification  of  Tlieorem  7  is  almost  identical  to  that  of 

Tlieorein  8,  wc  will  focus  only  on  the  first.  The  validity  of  Theorem  7 

is  suggested  by  noting  that  a  zero  (pole)  at  z=Zp  and  a  zero  (pole)  at 
★ 

z=1/Zq  contribute  the  same  magnitude  to  the  Fourier  transform.  Therefore, 
with  magnitude  information  alone,  tliere  is  an  inherent  ambiguity  in  the 
specification  of  the  sequence  in  that  a  zero  (pole)  of  the  original 
sequence  can  potentially  only  be  identified  from  the  magnitude  as  either 

■k 

a  zero  (pole)  at  or  at  z=1/Zq.  In  Theorem  7  this  ambiguity  is 

resolved  by  restricting  the  poles  to  1 ic  inside  the  unit  circle  and  the 
zeroes  to  lie  outside  while  in  Theorem  8,  the  condition  is  the  reverse. 

The  additional  condition  that  there  are  no  conjugate  reciprocal  pole-zero 
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pairs  eliminates  the  possibility  of  all-|)ass  terms  which  contribute  only 
to  the  phase  and  not  the  magnitude. 

io  more  formally  demonstrate  Theorem  7,  we  define  the  set  of 
sequences  to  consist  of  all  sequences  x[n]  which  satisfy  the  conditions 
of  Theorem  7  so  that  their  c-transform  are  restricted  to  be  of  the  form 

"l 

n  (l-a  cj 
k=l  ^ 

X(zJ  =  Cz"o  -  (17) 

n  (i-b,z  ) 
k=l 

■k 

with  for  any  k  and  f ,  and  where  C  is  a  real  constant,  n^  is  an 

integer,  and  |aj^|<l,  |bj^|<l  for  all  k.  Now  consider  two  sequences 
x[n]  and  y[n],  both  in  the  set  C^,  with  z-transforms  X(z)  and  Y(z) 
respectively.  We  wish  to  show  that  if  X(z)  and  Y(z)  both  have  the  same 
magnitude  on  the  unit  circle  then  x[n]=+y[n+m]  for  some  integer  m. 

Consider  G(z)  defined  as  the  ratio  X(z)/Y(z).  Since  X(z)  and  Y(z)  both 

have  the  same  magnitude  on  the  unit  circle,  G(z)  must  be  entirely  all-pass 

with  unity  magnitude,  i.e.,  for  a  zero  at  z=-q  there  must  be  a  pole  at 
★ 

z=1/Zq  and  vice  versa.  Therefore,  r,i.z)  consists  only  of  poles  and/or 

zeroes  at  z=0  or  at  z=°°  and  conjugate  reciprocal  pole-zero  pairs. 

Because  of  the  conditions  in  Theorem  7,  this  in  turn  requires  that  for 

any  zero  (or  pole)  of  X(z)  at  z=Zq  there  must  be  a  zero  (or  pole)  of 
★ 

Y(z)  at  z=l/Zp  which,  for  z^  /  0  or  violates  the  conditions  in 
Theorem  7  since  one  will  always  be  inside  and  the  other  outside  the 
unit  circle.  Thus,  G(z)  must  be  of  the  form 


or,  equivalently, 


for  some  integer  m. 


GCz)  =  +  z 


x[n]=+  y[n+m] 


A  conceptual  algorithm  similar  to  that  considered  in  Section  II 
can  be  developed  for  reconstruction  of  a  sequence  to  within  an  algebraic 
sign  and  a  delay  from  the  magnitude  of  its  Fourier  transform  under  the 
conditions  of  Theorem  7  or  Theorem  8.  We  outline  the  procedure  below 
under  the  conditions  of  Theorem  7.  It  is  easily  modified  for  the 
conditions  of  Theorem  8. 

Let  lx(w)|  denote  the  specified  magnitude  function.  Using  the 

Hilbert  transform,  a  minimum  phase  sequence  can  be  specified 

which  has  the  same  magnitude  function.  The  z-transform  X  .  (z)  of 

min 


X  ,  fn]  is  given  by 
min  ^  ^ 


X  .  (z)  =  C 

min  '  ' 


*  -1 

n  (l-a  z 
k=l 


(i-b.z‘^) 

k=l 


i.e.,  it  has  the  same  poles  as  X(z)  and  the  zeroes  are  reflected  inside 

the  unit  circle.  Since  the  conditions  of  Theorem  7  exclude  the 

* 

possibility  of  pole-zero  cancellation,  the  coefficients  a^,  and  thus 

a,  in  eq.  (17),  can  be  obtained  from  the  zeroes  of  X  .  (z)  and  the 

coefficients  b,  in  eq.  (17)  can  be  obtained  from  the  poles  of  X  .  (z) . 

k  ^  ^  min 
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V.  NUMIiRICAL  ALGORITHMS  FOR  RECONSTRUCTION  FROM  SAMPLES  OF  A 

PHASL  FUNCTION 

In  Section  II,  we  presented  two  sets  of  conditions,  embodied  in 
Theorems  1  and  2,  under  which  a  sequence  is  uniquely  specified  to 
within  a  positive  scale  factor  by  the  phase  of  its  Fourier  transform. 

In  this  section,  we  describe  two  numerical  algorithms  which  can  be  used 
to  reconstruct  a  sequence  satisfying  the  requirements  of  Theorem  1  from 
samples  of  its  phase  function  when  the  location  of  the  first  non-zero 
point  of  x[n]  is  known.  Although  these  algorithms  will  only  be  discussed 
in  terms  of  reconstructing  sequences  satisfying  the  conditions  of  Theorem  1, 
the  reconstruction  of  sequences  meeting  the  requirements  of  Theorem  2  may 
be  accomplished  by  simply  reconstructing  the  finite  length  sequence  x[n] 
defined  in  eq.  (5)  using  the  negative  of  the  specified  phase  samples  and 
then  computing  the  convolutional  inverse  sequence. 

The  first  algorithm  presented  below  is  an  iterative  technique  in 
which  the  estimate  of  x[n]  is  improved  in  each  iteration.  This  algorithm 
has  been  developed  by  slightly  modifying  an  iterative  technique  proposed 
by  Quatieri  (6)  for  reconstructing  a  minimum  phase  or  maximum  phase 
sequence  from  its  phase  function.  The  second  algorithm  is  a  closed  form 
solution  which  is  obtained  by  solving  a  set  of  linear  equations.  Under 
the  conditions  specified  in  Theorem  1,  this  algorithm  provides  the  desired 
sequence  x[n]  to  within  a  scale  factor  when  the  location  of  the  first 
non-zero  point  of  x[n]  is  known. 

In  the  discussions  which  follow,  x[n]  is  used  to  denote  a  sequence 
which  satisfies  the  conditions  of  Theorem  1  and  is  zero  outside  the 
interval  0<n<N-l  with  x[0]^0.  In  the  more  general  case  (see  Footnote  5), 
a  linear  phase  term  may  be  added  to  the  given  phase  to  accomplish  this. 
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A .  Iterative  M^’orithm 

The  M-poiiit  discrete  l-ourier  transform  (OFT)  of  x[n]  will  be  denoted 


as 

X(kJ  =  |X(k)  le-’^x-''^  (21) 

where  it  assumed  that  M>2N.  Then,  an  iterative  technique  to  reconstruct 
the  sequence  x(nj  from  the  M  samples  of  its  phase,  k=0, 1 , . .  .  ,M-1 ,  is 

as  follows: 


Step  1: 

We  begin  with  |x^(k)|,  an  initial  guess  of  the  unknown  OFT  magnitude 
and  form  the  first  estimate,  Xj(k),  of  X(k)  using  the  specified  phase 
function,  i.e, , 

X^(k)  =  lx^(k)|ej®x^''^  (22) 

Computing  the  inverse  DFT  of  Xj (k)  provides  the  first  estimate,  Xj[n], 
of  x[n].  Since  an  M-point  OFT  is  used,  Xj(n]  is  an  M-point  sequence 
which  is,  in  general,  non-zero  for  N<n<M-l. 


Step  2: 

From  Xj[ii],  another  sequence,  y^ln],  is  defined  by 


[nj 


x^[n] 

0 


0<n<N-l 

(23) 

N<n<M-l 
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Step  3: 


The  magnitude  |Y^(k)l  of  the  M-point  1>FT  of  y^tn]  is  then  considered 
as  a  new  estimate  of  |X(k}|  and  a  new  estimate  of  X(k)  is  formed  by 

X,(k)  =  |Yj(k)|e-'°x^''^  (24) 

From  this,  a  new  estimate  X2[nJ  is  obtained  from  the  inverse  DFT  of 
X^Ck).  Repetitive  application  of  steps  two  and  three  defines  the 
iteration. 

In  this  iterative  procedure,  the  total  squared  error  between  x[n] 

and  its  estimate  is  non-increasing  with  each  iteration.  To  see  this, 

let  x^[n]  denote  the  estimate  after  the  rth  iteration  and  define  the 

error,  E  ,  as 
r 

M-1 


|x[n]-x^[n] j 


(25) 


From  Parsevals's  theorem. 


1 


1 

M 


M-1 

Z 

k=0 


M-1 

Z 

k=0 


lx(k)-lx^(k)lej®x^’^^i 
|x(k)-lY^  (k)|e-'®x^^^ 

X  -  1 


(26) 


Since  Q^(k)  is  the  phase  of  X(k),  then 


lX(k)-lY^_^Ck)|e^®x^’'^ 


2  2 
<lX(k)-Y^_j(k)l 


(27) 


with  equality  if  and  only  if  0^(k)=0^  j^(k)  where  0^  ^(k)  is  the  phase 


17 


of  jllv).  fombininy  oq.  (.20)  and  (27)  gives 

M- 1  2 

>><M  "  lX(k)-Y  (k)l 

k  =  0 


ur,  with  I’arscval '  s  theorem; 


K  < 
r  - 


M-1 

Z 

n=0 


(28) 


Since 


V  , [nl  =  X  , In]  for  0  <n  <N-1 
'  r-1  ^  r-  1  ^ 


^  [n]  =  x[n]=0  for  N  <n  <M-1 


then, 


M-1 

Z 

n=0 


- 

1 

M-1 

- 

x[n]-y^_^[n] 

<  Z 
n=0 

xln]'X^-l[n] 

=  E 


r-1 


(29) 


with  equality  if  and  only  if  y^  ^[n]  =  x^  j[n]o  Therefore, 


E  <  E  , 
r  -  r-1 


(30) 


Although  eq.  (30)  is  not  sufficient  to  guarantee  the  existence 

/V  A 

of  a  converging  solution,  if  x[n]  is  a  converging  solution  then  x[n]= 

Ax[n]  with  A  >  0.  To  show  this,  we  note  that  the  equality  in  equation  (30) 
holds  if  and  only  if  0.^  j(k)=0^(k)  from  equation  (27)  and  if  and  only  if 
yj,_jfn]  =  x^  j(nl  from  equation  (29).  rhereforc.  if  the  algorithm  converges,  then 
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ii 


the  convergent  solution  x[nl  is  zero  outside  the  interval  t)<n<N-l  and 

n 

has  the  same  pl\ase  samples  as  xjuj  at  uj  =  k  tor  k=0, 1 ,  .  .  .  ,M- 1 .  Since 
M>2.'v,  there  exist  at  least  N-1  distinct  values  of  ic  in  the  interval 

at  whicli  equals  O^CcjJ  .  Then  from  Theorem  5,  x[n]=Ax[n]  for 

some  positive  constant  A. 

It  is  not  yet  known  if  this  iterative  procedure  always  leads  to  a 
converging  solutio)^.  However,  in  all  the  examples  that  we  have  considered 
so  far,  we  have  empirically  observed  that  the  algorithm  converges  to  the 
correct  solution  when  M>2N  and  that  the  number  of  iterations  required  to 
achieve  a  small  total  squared  error  is,  in  general,  quite  large.  We  have 
also  observed  tliat  increasing  M  may  increase  the  rate  of  convergence  of 
the  algoritlim,  but  such  an  increase  obviously  results  in  an  increase  in 
the  number  of  computations  required  for  each  iteration. 

Two  examples  of  the  iterative  procedure  applied  to  a  mixed  phase 
sequence,  x[n],  of  length  8  are  shown  in  Table  1.  In  the  first  example 
an  FFT  of  length  16  was  used.  In  the  second  example,  the  FFT  length  was 
extended  to  128  points.  In  both  cases,  the  initial  guess  of  the  unknown 
magnitude  was  chosen  to  be  a  constant,  and  the  scaling  factor  A  was 
chosen  so  that  the  resulting  sequences  have  the  same  value  at  the  origin 
as  x[nj.  The  results  after  10,  100,  500,  and  1000  iterations  are  presented 
along  with  the  values  of  the  total  squared  error. 

B.  Closed  ['orm  Solution 

4 

A  closed  form  solution  for  reconstructing  a  sequence  x[n]  from 
samples  of  its  phase,  G^(w),  follows  from  the  definition  of  G^(w) .  With 
defined  so  that  -7T<  0j^C<a)  <  't,  we  have 

4 

A  closed  form  solution  similar  to  the  one  presented  in  this  section 
can  be  obtained  by  expressing  the  real  and  imaginary  parts  of  X(iij)  in 
terms  of  the  given  phase  and  then  relating  the  real  and  imaginary 
parts  through  the  discrete  Hilbert  transform  relations  for  causal  sequences. 
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TABLE  1: 

Iterative  Reconstruction  of  a  Sequence  from  its  Phase 
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N-1 

-L  X  [  11 )  s i  n  nuj 

tan  (.ilj 

N-1 

I  x[nlcos  nu) 
n=() 


I'oi'  ttie  case  in  whicli  0_^(u)J=+rr/2  so  that  tan  0^(cu)=+«>,  eq.  (51)  is 
equiw'alcnt  to; 

N-1 

Z  x[  nj  cos  ncj  =  0  (52) 

n=0 


Sampling  tan  at  N-1  distinct  frequencies  Wj  ,  .  .  .  ^  with 

0  <  w,  <  ir  for  k=  1 , 2 ,  .  .  .  ,N-1 ,  eq.  (.51)  and  (52)  can  be  written  as: 

K 


N-1 

Z  X  [n] 
n=l 


cos  nwp,  tan  *  sin  nw^^ 


=  -x[0]  tan  if  0^ (Wj. ) 


(55)a 


N-1 

Z  x[nl  cos  nw  =  -x[0]  if  Q  (w, )=+Z_ 

,  K  X  K  -  T 


for  k=l , 2, . . .N- 1 .  tiquation  (35)  represents  (N-1)  linear  equations 
in  the  (N-l)  unknowns  of  x[nj  ard  can  be  expressed  in  matrix  form  as; 

S  X  =  x[0]b  (34) 


where  x  represents  the  vector  of  elements  of  x[n]  excluding  the  first 
element  of  tlie  sequence,  x[0] .  Any  solution  s[n],  to  eq.  (34)  has  the 
property  that  s[n]  is  zero  outside  the  interval  0<n<N-l  and  has  the  same 
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tanjiont  of  the  phase  as  x[nl  for  N-1  distinct  frequencies  between  zero 
and  71.  ihus,  from  i'heorem  5,  we  conclude  that  :  [nJ=Ax[nJ  for  some 
real  constant  A.  Iherefore,  there  is  only  one  jndejjendent  solution  to 
eq.  [.i4J.  This  implies  that  S  the  inverse  of  the  matrix  S,  exists 
and  that  x[n]  is  given  by 

X  =  AS’^h  i?,S) 

for  some  real  constant  A.  To  specify  x[nj  to  within  a  positive  scale 
factor,  we  first  assume  tliat  A>0  and  determine  the  phase  of  ^  in  eq.  (.^5). 
If  the  resulting  phase  does  not  differ  from  0^(u))  ,  then  A  is  positive. 
Otherwise,  A  is  negative. 

In  reconstructing  x[nj  from  using  eq.  (3S)  ,  it  should  be 

noted  that  we  have  some  control  over  the  matrix  S.  Since  the  elements 
of  the  matrix  S  are  functions  of  the  samples  of  G^(w) ,  S  can  be  changed 
by  choosing  a  different  set  of  frequency  samples.  This  control  over  S 
may  bo  useful  in  avoiding  those  frequency  samples  for  which  tan  is 

very  large  or  in  avoiding  potential  numerical  instabilities  in  computing 

.  .  .  0-1 

the  inverse  matrix  S 

Compared  with  the  iterative  algorithm,  the  closed  form  solution 
presented  above  has  the  advantage  that  the  desired  sequence  is  guaranteed 
to  be  the  solution  to  eq.  (.34)  and,  in  addition,  no  iterations  are 
required  in  order  to  reach  a  solution.  On  the  other  hand,  equation 
(54)  requires  computing  the  inverse  of  an  (N-l)x(N-l)  matrix  which  may 
lead  to  numerical  problems  and  severe  round-off  errors,  particularly  as 

I 

N  becomes  large. 

The  algoritlim  discussed  above  lias  been  applied  to  a  variety  of 
different  examples.  Consistent  with  our  theoretical  results,  in  all 
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examples  we  have  considered,  the  desired  solutions  have  been  obtained. 
Specifically,  for  the  sequence  shown  in  Table  1,  the  closed  form  solution 
was  used  to  reconstruct  the  sequence  from  its  phase.  The  phase  samples 
used  to  define  the  matrix  S  were  first  cliosen  to  be  equally  spaced 
between  zero  and  tt,  and  then  were  randomly  selected  within  the  interval 
zero  to  II.  Within  the  limits  of  finite  jirecision  arithmetic,  in  both 
cases  the  sequence  was  reconstructed  exactly  when  the  scaling  factor  was 
chosen  so  that  the  solutions  obtained  had  the  correct  value  at  the 
origin. 

VI.  EXTENSION  TO  MULTI-DIMTNSIONAL  SEQUENCES 

In  this  section,  we  extend  the  results  of  Sections  II,  III,  IV  and 
V  to  the  case  of  multi-dimensional  sequences.  This  extension  is  achieved 
by  mapping  a  multi-dimensional  sequence  into  a  1-D  (one  dimensional) 
sequence  and  then  applying  the  results  for  1-D  sequences.  Since  the 
extension  of  the  2-D  case  to  sequences  of  higher  dimension  is  straight¬ 
forward,  our  discussions  in  the  section  will  concentrate  on  the  2-D 
case.  Again,  we  consider  only  sequences  which  are  real,  have  rational  z- 
transforms,  and  have  Fourier  transforms  that  converge. 

Let  x[n^,n2]  represent  a  2-D  sequence  which  has  a  rational  z- 
transform  with  X(a)^,W2)  given  by 

A(a)  ) 

=  ___ 

where  ACco^,^^)  is  a  2-D  polynomial  of  degree  in  exp[j(i)j]  and  in 

exp[jw^]  and  where  is  a  2-D  polynomial  of  degree  in  expfjWj] 

and  in  exp[jcj2].  Suppose  we  form  a  1-D  sequence  x^[n]  or  X2[n]  by 

^  ^ 

x^[nj  =  E  x[m,n-Nm]  (37)a 
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> 


or  2  x[n-Mni,mJ  (37)b 

ni= 

vvhcre  M  >  maxCNj.M^J  and  N  >  max(M,,N,).  Then  it  can  be  shown  (7,8) 
t!\at  the  transformation  in  equation  (37)  is  invertible  and  that 

X^(u)  = 


(38)a 

Wj^=wN,  cj^=cj 


and 


X,(,u))  =  X(w^,u,) 


W2=wM 


(38)b 


From  equation  (38),  it  is  clear  that  the  phase  of  x^ [n]  or  X2[n]  is 

/'  /N 

specified  by  the  phase  of  xlUj^.n^]  and  the  magnitude  of  x^[n]  or  ^2^^^ 

is  specified  by  the  magnitude  of  x[nj,n2].  Therefore,  all  the  theorems 

and  numerical  algorithms  developed  in  Secitons  II,  III,  IV  and  V  for  1-D 

sequences  may  be  extended  to  2-D  sequences  by  first  transforming  them 

into  1-D  sequences  using  equation  (37)  and  then  applying  the  1-D  results 

to  the  resulting  1-D  sequences.  Thus,  for  example.  Theorem  1  may  be 

extended  to  2-D  sequences  as  follows:  Let  x[n^,n2]  and  y[n^,n2]  be  two 

2-D  sequences  such  that  the  two  1-D  sequences  x[n]  and  y[n]  obtained 

from  transforming  x[n^,n2J  and  y[nj,n2]  using  equation  (37)  are  finite 

in  length  and  have  no  zeroes  in  reciprocal  pairs  or  on  the  unit  circle. 

If  G  ((j  ,w^)  =  0  (cj  ,w  )  for  all  to  and  to  ,  then  x[n,  ,n,]=Ay [n  ,n_]  for 
Xl^  yl*:  1  /  i4 

some  positive  constant  A.  If  tan  G^(t0j,u)2)  =  tan  0y(tOj^,tO2)  for  all 
(Oj  and  then  x[nj^,n2]=Ay[n^,n2]  for  some  real  constant  A. 

If  x[n^,n^J  is  a  2-D  sequence  with  finite  support,  then  the  transformation 
given  by  equation  (57)  can  be  reduced  to  a  simpler  form.  Specifically, 
let  x[nj,n2]  be  zero  outside  the  region  0<n^<N^-l  and  0<n.,<N2-l.  Then 
equation  (37)  can  be  rewritten  as 
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x^fiijN+n^l  =  x[nj,n,J  with  N>Nj 


(39)a 


or  x,[ii^+n^M]  =  x|n^,n,,j  with  M>N2  (39Jb 

Clearly,  the  transformation  in  equation  (.39)  is  invertible  and  it  can  be 
easily  shown  (8)  that  Xj(tjj)  and  X^Cai)  are  given  by  equation  (38). 

As  an  illustration  of  the  results  of  this  section,  a  2-D  sequence 

representing  the  intensity  of  an  image,  x[n^,n2],  was  created  which  is 
zero  outside  the  region  0<nj^<ll  and  0<n2<ll.  From  the  phase  of  Xj^(w) 
defined  by  equation  (38)a  with  .N'=12,  the  closed  form  solution  was  used 
to  reconstruct  x[n^,n,].  With  the  scale  factor  chosen  so  that  the 
reconstructed  image  had  the  same  value  at  the  origin  as  x[nj,n2],  the 
result,  shown  in  Figure  1,  is  indistinguishable  from  the  original.  For 
illustration,  the  image  shown  has  been  enlarged  by  means  of  a  zero-order 
hold. 

Finally,  we  present  another  example  which  illustrates  a  potential 
application  of  the  results  of  this  paper  to  blind  deconvolution.  In 
ireage  processing,  a  model  which  frequently  arises  to  describe  image  blur 
in  the  absence  of  noise  is: 

y[n^,n2]  =  x[nj,n2]  *  htn^,n2]  (40) 

where  x[nj^,n2J  is  the  original  unblurred  image  and  h[n^,n2]  is  a  (real) 
symmetric  2-D  sequence,  i.e.,  h[nj,n2]  =  h[-nj,-n2].  This  model  may  be 
used,  under  the  appropriate  band-limited  assumptions,  to  account  for 
uniform  motion  blur,  severely  defocussed  lenses  with  circular  aperature 
stops,  and  long-term  exposure  of  atmospheric  turbulence  [9].  In  the 
context  of  blind  deconvolution,  the  desired  objective  is  to  recover 
x[nj,n2]  in  eq.  (40)  from  the  blurred  image  y[n^,n2]  with  no  detailed 
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knowledge  of  However,  if  hln^.n^l  is  known  to  be  symmetric, 

one  possible  approach  to  the  deconvolution  of  eq,  (40)  is  suggested  by 
the  results  of  this  section  and  section  V.  S])ec i f i cal ly ,  let  x^n^.n^] 
denote  an  original  Linl)lurred  image  and  liln^.n,)  a  symmetric  point  spead 
function.  Witli  >'[n^,n,]  rejiresent  ing  tlie  lilurred  image  as  expressed  by 
equation  (40) ,  define  the  coordinate  system  for  yln^^.n^]  in  such  a  way 
that  y[n^,n.,J  is  zero  outside  the  region  ()vnj"  .\|-l  and  i)<n.,<N.,-l .  Since 
h[nj,n,J  is  symmetric,  its  phase  is  zero  or  z  for  all  i.jj  and  . 
Therefore, 

tan 

for  all  ojj^  and  oj., .  Suppose  there  exists  a  sequence  Xj[nl  or  x^fn] 
defined  by  equation  (39)  which  satisfies  the  conditions  of  Theorem  3 
and  that  the  location  of  the  first  non-zero  point  of  this  sequence  is 
known.  To  simplify  the  following  discussion,  we  denote  this  sequence 

/S  /N 

simply  by  x[n]  and  define  the  sequences  y[n]and  h[n]  by  j^erforming  the 
same  transformation  on  y[n^,n2]  and  h[nj,n2]  that  is  used  to  obtain 
x[n] .  From  the  results  of  (8) , 

/N  /S  /S 

y[n]  =  x[n]  *  h[n]  (42) 

and  therefore,  it  follows  from  equations  (38)  and  (41)  that 

tan  0/v(a))  =  tan  0.^(0))  (43) 

y  X 

A, 

for  all  u.  Since  the  location  of  the  first  non-zero  point  of  x[n]  is 
assumed  to  be  known,  by  adding  the  appropriate  linear  phase  to  0js(w), 


2b 


Image  reconstruction  from  phase  information. 
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Fig.l.  Image  reconstruction  from  phase  information. 


the  corresponding  sequence,  x[nj,  satisfies  the  conditions  of  Theorem  5, 
i.e.,  x[nj  satisfies  the  conditions  of  Theorem  3  with  x[0J^0.  Consequently, 
the  closed  form  solution  given  by  equation  (55)  specifies  x[n]  to 
within  a  scale  factor  from  samples  of  its  phase  or,  equivalently,  from 
samples  of  the  phase  of  yln^.n^]-  Therefore,  since  equation  (39)  represents 
an  invertible  transformation,  x[nj,n2J  can  be  recovered  from  x[n]  and 
knowledge  of  the  location  of  the  first  non-zero  point  of  x[n]. 

Since  it  is  unreasonable  to  assume  that  the  location  of  the  first 
non-zero  point  of  x[n]  is  always  known,  we  outline  one  possible  precedure 
for  finding  this  point  from  y[n]  which  is  given  by  equation  (42) . 
Specifically,  with  y[n]  equal  to  zero  outside  the  interval  0<n<N-l, 
suppose  we  assume  that  x[0]^0  and  proceed  to  solve  equation  (34)  for  the 

A 

vector  _x  using  (N-1)  samples  of  0^(a)) .  Then,  it  is  easy  to  show  that 
the  (N-1)  X  (N-1)  matrix  S  in  equation  (34)  is  singular  if  x[0]^0. 

More  importantly,  however,  it  can  be  shown  that  n^,  the  location  of  the 
first  non-zero  point  of  x[n],  is  related  to  the  rank  of  S  by: 

n^  =  (N-1)  -  rank[S]  (44) 

Therefore,  suppose  we  are  given  the  2-D  sequence  y[n^,n2]  in  equation 
(40)  with  y[n^,n2]  taken  to  be  zero  outside  the  region  0<nj<N^  and 

A  A 

0<n2<n2-l.  Then,  if  x^[n]  or  x^fn]  satisfy  the  conditions  of  Theorem  3, 
equation  (44)  may  be  used  to  determine  the  location  of  the  first  non- 

A  A 

zero  point  of  x^ [n]  or  X2[n].  With  this  information,  equation  (34)  may 

A  A 

be  solved  to  recover  x^[n]  or  x^in]  to  within  a  scale  factor  from  ©^(u)^,^^). 
Then,  performing  the  inverse  of  the  transformation  given  by  equation 
(39),  x[n^,n2J  can  be  reconstructed  to  within  a  scale  factor. 
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An  example  which  illustrates  the  results  of  this  section  to  blind 
deconvolution  is  shown  in  I'igure  2,  The  image  in  Figure  2a  has  a 
support  of  lb  pels  x  lb  pels  whicli,  for  illustration,  has  been  enlarged 
using  a  zero-order  hold  and  represents  the  original  image,  x[n^,n^l  , 
w'hicti  has  been  blurred  by  a  symmetrical  2-1)  sequence  h[iij,n_,J  as  in 
eq.  (40).  After  determining  the  suj)port  of  h[npn2J  from  the  blurred 
image  using  a  variation  of  the  procedure  outlined  above  for  determining 
the  location  of  the  first  iion-zero  point  of  x^|n]  or  x.,fn],  the  image 
shown  in  Figure  2b  was  obtained  from  the  closed  form  solution.  Within 
the  finite  precision  of  digital  computation,  the  reconstruction  is  exact. 

Although  the  iterative  technique  may  also  be  used  in  blind  deconvolution, 
it  has  an  important  limitation  as  compared  to  the  closed  form  solution. 
Specifically,  the  iterative  algorithm  requires  that  the  phase  of  the 
blurred  image  be  identical  to  the  phase  of  the  original  image  as  opposed 
to  the  tangents  of  the  phase  being  equal  as  required  in  the  case  of  the 
closed  form  solution.  This,  in  turn,  restricts  the  blurring  function 
hfn^,n,J,  to  be  zero  phase. 

Although  we  have  presented  one  possible  approach  to  blind  deconvolution 
when  the  degrading  function  is  symmetric,  it  does  not  constitute  a 
general  solution  to  the  problem  of  blind  deconvolution  since,  in  practice, 
eq.  (40)  is  only  an  approximation.  More  generally,  a  blurred  image  is 
given  by 

=  |x[n^,n2]  *  h[nj .n,] |  w[n^ .n^J  +  v[n^,n2]  (46) 

where  v[n^,n2]  is  additive  noise  and  w[nj,n2]  is  a  window  function.  The 
sensitivity  of  our  proposed  deconvolution  technique  to  the  approximation 
in  eq.  (40)  as  well  as  possible  extensions  to  the  more  general  case  in 
eq.  (46)  are  currently  under  investigation. 
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Fig,2(.a-b),  Blind  deconvolution  of  images:  la)  Image  blurred  by  symme 
point  spread  function,  (b)  Result  of  blind  deconvolution  of  image  (a) , 
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